** Clearing Stata memory
capture log close
clear all
set more off, perm
set seed 1234

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////// Table O.43: Alternative Residuals: Separate regressions by subject and gender  ///////////////
//////////////////////////// Control for all P1 subject scores, Log (Annual Wages) Seven to 12 Years After Admission Exam /
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

** Opening Phase 2 norm_scores dataset 
use "Work Data/Gender_Phase2_long.dta",clear

*** Creating variables
encode subject, gen (sub)
tab subject, gen (d_sub)
label var sub "Subject"

** Subject dummies
rename d_sub1 Biology
rename d_sub2 Chemistry
rename d_sub3 Geography
rename d_sub4 History
rename d_sub5 Language
rename d_sub6 Mathematics
rename d_sub7 Physics
rename d_sub8 Portuguese
* Labels
label var Biology "Biology"
label var Chemistry "Chemistry"
label var Geography "Geography"
label var History "History"
label var Math "Mathematics"
label var Physics "Physics"
label var Portuguese "Portuguese"
label var Language "Foreign Language"

* Days of admission exam
gen exam_day=1 if Portuguese==1 | Biology==1 //Day 1: Portuguese and Biology
replace exam_day=2  if Chemistry==1 | History==1 //Day 2: Chemistry and History
replace exam_day=3  if Physics==1 | Geography==1 //Day 3: Physics and Geography
replace exam_day=4  if Mathematics==1 | Language==1 //Day 4: Mathematics and English
tab exam_day
label var exam_day "Days of admission exam"

* # of priority discipline per applicant, per day of the admission exam
levelsof exam_day, local(levels) 
foreach i of local levels {
bys inscri2: egen n_prior_day`i'=sum(priority) if exam_day==`i'
bys inscri2: egen n_priority_day`i'=max(n_prior_day`i')
drop n_prior_day`i'
sum n_priority_day`i'
label var n_priority_day`i' "\# priority subjects in exam day `i'"
}
tab career_choice if n_priority_day1==2 // only 1 career has 2 priority subjects in the same day of the admission exam (Phonology)

** Create other priority variables:

* Only a given subject is priority:
gen only_priority=0
levelsof exam_day, local(levels) 
foreach i of local levels {
replace only_priority=1 if priority==1 & exam_day==`i' & n_priority_day`i'==1
}

* Only the other subject is priority:
gen other_priority=0
levelsof exam_day, local(levels) 
foreach i of local levels {
replace other_priority=1 if priority==0 & exam_day==`i' & n_priority_day`i'==1
}

* Both priority
gen both_priority=0
levelsof exam_day, local(levels) 
foreach i of local levels {
replace both_priority=1 if priority==1 & exam_day==`i' & n_priority_day`i'==2
}

** P1 scores: P1 normalized subject-specific scores
forvalues i=2(1)4 {
gen norm_p1score`i'=norm_p1score^`i'
sum norm_p1score`i'
}

*********************************************************************************
****************   Relative performances ****************************************
*********************************************************************************

** ENEM

foreach v in norm_enem_w {
bys year female: egen `v'_ave_g=mean(`v')
gen `v'_g=`v'-`v'_ave_g
bys year female: sum `v'_g
}
drop norm_enem_w_ave_g

** Phase 1 scores

foreach v in norm_p1score {
tab year, sum(`v')
bys year subject female: egen gs_`v'_ave=mean(`v')
gen gs_`v'=`v'-gs_`v'_ave
bys year female subject: sum gs_`v'
drop gs_`v'_ave
}

*********************************************************************************
**************** Main sample ****************************************************
*********************************************************************************

* 1) Only years before the affirmative action took place
drop if aa_year==1
drop if year==2000
tab year

 ******************************** SAVE DATA SET TO BOOTSTRAP ********************************
 
keep norm_score gs_norm_p1score norm_enem_w_g subject only_priority other_priority both_priority sex female norm_enem_w career_choice year inscri2 priority

save "Work Data/Temp_Bootstrap_AddTab.dta", replace

*********************************************************************************
****************   Regressions by exam day **************************************
*********************************************************************************

sort inscri2
order subject gs_norm_p1score 

local j=1	
levelsof subject, local(levels) 
foreach i of local levels {

	by inscri2: egen norm_p1score_`i' = mean(gs_norm_p1score) if subject=="`i'" // relative performance in phase 1, similar to main specification
	by inscri2: egen aux_`i' = max(norm_p1score_`i') 
	replace norm_p1score_`i' = aux_`i' 
	drop aux_`i' 
	
	if `j' ==1 loc title= "Bio"
	if `j' ==2 loc title= "Chem"
	if `j' ==3 loc title= "Geog"
	if `j' ==4 loc title= "Hist"
	if `j' ==5 loc title= "EngFr"
	if `j' ==6 loc title= "Math"
	if `j' ==7 loc title= "Phys"
	if `j' ==8 loc title= "Port"
		
	label var norm_p1score_`i' "P1 scores `title' "

	local ++j
}

global norm_p1score_all norm_p1score_biol norm_p1score_chem norm_p1score_geog norm_p1score_hist norm_p1score_math norm_p1score_phy norm_p1score_port

estimates clear

global outcome_list = "norm_score"

local count_outcome = 1
foreach outcome in $outcome_list {

if `count_outcome' ==1 local outcome_name = "scores"
if `count_outcome' ==1 local outcome_title = "Scores"

local count_enem=1

foreach j in norm_enem_w_g {

	levelsof subject, local(levels) 
	foreach i of local levels {
	di "`i'" 
	di "`j'"
	
	forvalues sex =1/2{ // 1= male, 2=female
	
		if "`i'"!="lang"{

		reg `outcome' only_priority other_priority both_priority  `j'  $norm_p1score_all if subject=="`i'" & sex==`sex', noomitted
		* Saving residuals
		 predict resall_`sex'_`i', residuals

			}
		} 
	}

		local ++count_enem
}
local ++count_outcome
}


**** Constructing measures of the residuals

drop if subject=="lang" 

gen resall=.
foreach x in 1 2 {
levelsof subject, local(levels) 
foreach s of local levels {
replace resall=resall_`x'_`s' if subject=="`s'" & sex==`x'
}	
}
sum resall

foreach v of varlist resall {
bys inscri2: egen `v'_prio_ave=mean(`v') if priority==1
bys inscri2: egen `v'_prio_average=min(`v'_prio_ave)
mdesc `v'_prio_average
bys inscri2: egen `v'_nonprio_ave=mean(`v') if priority==0
bys inscri2: egen `v'_nonprio_average=min(`v'_nonprio_ave)
mdesc `v'_nonprio_average
gen `v'_diff=`v'_prio_average-`v'_nonprio_average
mdesc `v'_diff
ttest `v'_diff, by(female)
}

collapse (mean) resall_diff female norm_enem_w career_choice year, by(inscri2)

 ***************************************************
 ****************** Wages RAIS *********************
 ***************************************************

count
merge 1:1 inscri2 using "Work Data/RAIS_cleaned.dta"
drop if _merge==2
tab _merge 

forvalues i=7(1)12 {
gen yearafter`i'=year+`i'
tab yearafter`i', mi
gen annual_wage_after`i'=.
levelsof yearafter`i', local(levels) 
foreach l of local levels {
replace annual_wage_after`i'=mwagetot`l' if yearafter`i'==`l'
sum  annual_wage_after`i' 
}
}

**** Average and maximum wage

foreach x in annual_wage {
    
* Average
egen avg_`x'_712=rowmean(`x'_after7 `x'_after8 `x'_after9 `x'_after10 `x'_after11 `x'_after12) // 7-12 years after admission exam
sum avg_`x'*
sum avg_`x'* if _merge==1

* Maximum
egen max_`x'_712=rowmax(`x'_after7 `x'_after8 `x'_after9 `x'_after10 `x'_after11 `x'_after12) // 7-12 years after admission exam
sum max_`x'*
sum max_`x'* if _merge==1
}

* Log variables
foreach v of varlist annual_wage_after*  max* avg*  {
 gen l_`v'=log(`v')
 sum l_`v'
}

foreach v of varlist res*diff {
sum `v'
gen mean_`v'=r(mean)
gen sd_`v'=r(sd)
gen norm_`v'=(`v'-mean_`v')/sd_`v'
sum  norm_`v'
}

******************************** SAVE WAGES DATA TO MERGE WITH DATA SET TO BOOTSTRAP ********************************
preserve
keep l_*_annual_wage_712 inscri2 norm_*diff career_choice norm_enem_w female year
renvarlab norm_*diff career_choice norm_enem_w female, postfix(_orig)  
duplicates report
save "Work Data/Temp_Wages_AddTable_BS.dta", replace
restore

foreach m of varlist norm_*diff {

** # 5 - Control for the residuals 
drop if `m'==.

**********  Annual Wages and Participation Formal Labor Market **********

foreach v of varlist l_*_annual_wage_712 {

preserve

capture erase "Work Data/BS_addtab_allp1_`m'_`v'.dta"

** Do not control for program FE

reghdfe `v' female, vce(robust) absorb(year)
foreach x in _cons female {
gen b_`x'1=_b[`x']
gen se_`x'1=_se[`x']	
}

reghdfe `v' female `m', vce(robust) absorb(year)
foreach x in _cons female `m' {
gen b_`x'2=_b[`x']
gen se_`x'2=_se[`x']	
}

reghdfe `v' female norm_enem_w, vce(robust) absorb(year)
foreach x in _cons female norm_enem_w  {
gen b_`x'3=_b[`x']
gen se_`x'3=_se[`x']	
}

reghdfe `v' female `m' norm_enem_w, vce(robust) absorb(year)
foreach x in _cons female `m' norm_enem_w  {
gen b_`x'4=_b[`x']
gen se_`x'4=_se[`x']	
}

** Control for program FE

reghdfe `v' female, vce(robust) absorb(career_choice year) 
foreach x in _cons female {
gen b_`x'5=_b[`x']
gen se_`x'5=_se[`x']	
}

reghdfe `v' female `m', vce(robust) absorb(career_choice year)
foreach x in _cons female `m' {
gen b_`x'6=_b[`x']
gen se_`x'6=_se[`x']	
}

reghdfe `v' female norm_enem_w, vce(robust) absorb(career_choice year)
foreach x in _cons female norm_enem_w  {
gen b_`x'7=_b[`x']
gen se_`x'7=_se[`x']	
}

reghdfe `v' female `m' norm_enem_w, vce(robust) absorb(career_choice year)
foreach x in _cons female `m' norm_enem_w  {
gen b_`x'8=_b[`x']
gen se_`x'8=_se[`x']	
}

gen boot_run=0
keep boot_run b_* se_*
keep if _n==1  
count
sum

save "Work Data/BS_addtab_allp1_`m'_`v'.dta", replace

restore

}
}


********************************************************************************************************
**************************************  Step 2: Bootstrap   ********************************************
********************************************************************************************************

use "Work Data/Temp_Bootstrap_AddTab.dta", clear
merge m:1 inscri2 using "Work Data/Temp_Wages_AddTable_BS.dta"
drop _merge
count 

renvarlab norm_resall_diff_orig , postsub("_orig" "")

local boots = 999

quietly {
forvalues b=1(1)`boots' {

noisily display "Working on `b' out of `boots' at $S_TIME"

foreach m of varlist norm_*diff {
	
foreach v of varlist l_*_annual_wage_712 {

preserve

local seed=2023 + `b'
set seed `seed'

bsample,  cluster(inscri2) idcluster(newid)  

********** Estimating the residuals **********

global outcome_list = "norm_score"

local count_outcome = 1
foreach outcome in $outcome_list {
local count_enem=1

foreach j in norm_enem_w_g {

	levelsof subject, local(levels) 
	foreach i of local levels {
	di "`i'" 
	di "`j'"
	
	forvalues sex =1/2{ // 1= male, 2=female
	
		if "`i'"!="lang"{

			reg `outcome' only_priority other_priority both_priority  `j' gs_norm_p1score if subject=="`i'" & sex==`sex', noomitted
		* Saving residuals
		 predict resall_`sex'_`i', residuals

			}
		} 
	}

		local ++count_enem
}
local ++count_outcome
}

**** Constructing measures of the residuals

drop if subject=="lang" 

gen resall=.
foreach x in 1 2 {
levelsof subject, local(levels) 
foreach s of local levels {
replace resall=resall_`x'_`s' if subject=="`s'" & sex==`x'
}	
}
sum resall

foreach r of varlist resall {
bys newid: egen `r'_prio_ave=mean(`r') if priority==1
bys newid: egen `r'_prio_average=min(`r'_prio_ave)
mdesc `r'_prio_average
bys newid: egen `r'_nonprio_ave=mean(`r') if priority==0
bys newid: egen `r'_nonprio_average=min(`r'_nonprio_ave)
mdesc `r'_nonprio_average
gen `r'_diff=`r'_prio_average-`r'_nonprio_average
mdesc `r'_diff
ttest `r'_diff, by(female)
}

collapse (mean) res*diff female norm_enem_w year career_choice `v', by(newid)

* Normalizing residuals
foreach r of varlist res*diff {
sum `r'
gen mean_`r'=r(mean)
gen sd_`r'=r(sd)
gen norm_`r'=(`r'-mean_`r')/sd_`r'
sum  norm_`r'
}
	
drop if `m'==.

** Do not control for program FE

reghdfe `v' female, vce(robust) absorb(year)
foreach x in _cons female {
gen b_`x'1=_b[`x']
}

reghdfe `v' female `m', vce(robust) absorb(year)
foreach x in _cons female `m' {
gen b_`x'2=_b[`x']
}

reghdfe `v' female norm_enem_w, vce(robust) absorb(year)
foreach x in _cons female norm_enem_w  {
gen b_`x'3=_b[`x']
}

reghdfe `v' female `m' norm_enem_w, vce(robust) absorb(year)
foreach x in _cons female `m' norm_enem_w  {
gen b_`x'4=_b[`x']
}

** Control for program FE
reghdfe `v' female, vce(robust) absorb(career_choice year) 
foreach x in _cons female {
gen b_`x'5=_b[`x']
}

reghdfe `v' female `m', vce(robust) absorb(career_choice year)
foreach x in _cons female `m' {
gen b_`x'6=_b[`x']
}

reghdfe `v' female norm_enem_w, vce(robust) absorb(career_choice year)
foreach x in _cons female norm_enem_w  {
gen b_`x'7=_b[`x']
}

reghdfe `v' female `m' norm_enem_w, vce(robust) absorb(career_choice year)
foreach x in _cons female `m' norm_enem_w  {
gen b_`x'8=_b[`x']
}

gen boot_run=`b'
keep boot_run b_*
keep if _n==1  
summarize

append using "Work Data/BS_addtab_allp1_`m'_`v'.dta" 
save "Work Data/BS_addtab_allp1_`m'_`v'.dta", replace


restore

}


}
}
}


foreach m in norm_resall_diff {

********************************************************************************************************
*********************** Step 3: Generating the bootstrapped standard errors ****************************
********************************************************************************************************

use "Work Data/Temp_Wages_AddTable_BS.dta", clear

matrix drop _all

foreach y of varlist l_*_annual_wage_712 {
	
use  "Work Data/BS_addtab_allp1_`m'_`y'.dta", clear

** Calculating bootstrapped standard errors

foreach x in female norm_enem_w `m' _cons {
forvalues i = 1(1)8 {
capture egen BS_sd_`x'`i'=sd(b_`x'`i')  if boot_run~=0
}
}

** Making a matrix with bootstrapped standard errors

matrix define bsemat`y'1 = J(1,2,0)
matrix define bsemat`y'2 = J(1,3,0)
matrix define bsemat`y'3 = J(1,3,0)
matrix define bsemat`y'4 = J(1,4,0)
matrix define bsemat`y'5 = J(1,2,0)
matrix define bsemat`y'6 = J(1,3,0)
matrix define bsemat`y'7 = J(1,3,0)
matrix define bsemat`y'8 = J(1,4,0)

matrix colnames bsemat`y'1 = female _cons
matrix colnames bsemat`y'2 = female resid _cons
matrix colnames bsemat`y'3 = female norm_enem_w _cons
matrix colnames bsemat`y'4 = female resid norm_enem_w _cons
matrix colnames bsemat`y'5 = female _cons
matrix colnames bsemat`y'6 = female resid _cons
matrix colnames bsemat`y'7 = female norm_enem_w _cons
matrix colnames bsemat`y'8 = female resid norm_enem_w _cons

forvalues i = 1(1)8 {
sum BS_sd_female`i'
scalar BS_sd_female = r(mean)
matrix bsemat`y'`i'[1,1] = BS_sd_female
}

foreach i in 2 4 6 8 {
sum BS_sd_`m'`i'
scalar BS_sd_`m' = r(mean)
matrix bsemat`y'`i'[1,2] = BS_sd_`m'
}

foreach i in 3 7 {
sum BS_sd_norm_enem_w`i'
scalar BS_sd_norm_enem_w = r(mean)
matrix bsemat`y'`i'[1,2] = BS_sd_norm_enem_w
}

foreach i in 4 8 {
sum BS_sd_norm_enem_w`i'
scalar BS_sd_norm_enem_w = r(mean)
matrix bsemat`y'`i'[1,3] = BS_sd_norm_enem_w
}

foreach i in 1 5 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,2] = BS_sd_cons
}

foreach i in 2 6 3 7 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,3] = BS_sd_cons
}

foreach i in 4 8 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,4] = BS_sd_cons
}

forvalues i = 1(1)8 {
matrix list bsemat`y'`i'
}

}

********************************************************************************************************
*********************** Step 4: Computing p-values and exporting tables ********************************
********************************************************************************************************

use "Work Data/Temp_Wages_AddTable_BS.dta", clear
rename `m'_orig resid
rename career_choice_orig career_choice
rename norm_enem_w_orig norm_enem_w
rename female_orig female

label var norm_enem_w "Norm. ENEM scores"
label var female "Female"
label var resid "Relative priority performance"

**********  Annual Wages **********

estimates clear

drop if resid==.

foreach y of varlist l_*_annual_wage_712 {

** Do not control for program FE

reghdfe `y' female, vce(robust) absorb(year)
matrix define b`y'1 = e(b)
scalar DF`y'1=e(df_r)
estimates store `y'1
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'1
estadd matrix seboot =  bsemat`y'1

reghdfe `y' female resid, vce(robust) absorb(year)
matrix define b`y'2 = e(b)
scalar DF`y'2=e(df_r)
estimates store `y'2
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'2
estadd matrix seboot =  bsemat`y'2

reghdfe `y' female norm_enem_w, vce(robust) absorb(year)
matrix define b`y'3 = e(b)
scalar DF`y'3=e(df_r)
estimates store `y'3
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'3
estadd matrix seboot =  bsemat`y'3

reghdfe `y' female resid norm_enem_w, vce(robust) absorb(year)
matrix define b`y'4 = e(b)
scalar DF`y'4=e(df_r)
estimates store `y'4
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'4
estadd matrix seboot =  bsemat`y'4

** Control for program FE

reghdfe `y' female, vce(robust) absorb(career_choice year)
matrix define b`y'5 = e(b)
scalar DF`y'5=e(df_r)
estimates store `y'5
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'5
estadd matrix seboot =  bsemat`y'5

reghdfe `y' female resid, vce(robust) absorb(career_choice year)
matrix define b`y'6 = e(b)
scalar DF`y'6=e(df_r)
estimates store `y'6
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'6
estadd matrix seboot =  bsemat`y'6

reghdfe `y' female norm_enem_w, vce(robust) absorb(career_choice year)
matrix define b`y'7 = e(b)
scalar DF`y'7=e(df_r)
estimates store `y'7
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'7
estadd matrix seboot =  bsemat`y'7

reghdfe `y' female resid norm_enem_w, vce(robust) absorb(career_choice year)
matrix define b`y'8 = e(b)
scalar DF`y'8=e(df_r)
estimates store `y'8
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'8
estadd matrix seboot =  bsemat`y'8

estadd local program "Yes": `y'5  `y'6 `y'7  `y'8
estadd local program "No": `y'1  `y'2 `y'3  `y'4
estadd local year "Yes": `y'*

* Calculating p-values

foreach x in pval pvalue {	
	
matrix define b`x'`y'1 = J(1,2,0)
matrix define b`x'`y'2 = J(1,3,0)
matrix define b`x'`y'3 = J(1,3,0)
matrix define b`x'`y'4 = J(1,4,0)
matrix define b`x'`y'5 = J(1,2,0)
matrix define b`x'`y'6 = J(1,3,0)
matrix define b`x'`y'7 = J(1,3,0)
matrix define b`x'`y'8 = J(1,4,0)

matrix colnames b`x'`y'1 = female _cons
matrix colnames b`x'`y'2 = female resid _cons
matrix colnames b`x'`y'3 = female norm_enem_w _cons
matrix colnames b`x'`y'4 = female resid norm_enem_w _cons
matrix colnames b`x'`y'5 = female _cons
matrix colnames b`x'`y'6 = female resid _cons
matrix colnames b`x'`y'7 = female norm_enem_w _cons
matrix colnames b`x'`y'8 = female resid norm_enem_w _cons

}

foreach x in 1 5 {
forvalues j = 1/2 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

foreach x in 2 3 6 7 {
forvalues j = 1/3 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

foreach x in 4 8 {
forvalues j = 1/4 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

forvalues i = 1(1)8 {
matrix list b`y'`i'
matrix list bsemat`y'`i'
matrix list bpval`y'`i'
matrix list bpvalue`y'`i'
}

* Adding my p values to the estimation results.
forvalues i=1(1)8 {
estimates restore `y'`i'
estadd matrix pboot = bpval`y'`i'
}

}

* Average and maximum wages between 7-12 years after admission exam

esttab l_avg_annual_wage_712*  using "Output\Add_Table_AllP1_Annual_wages_`m'_ENEM_712_BS.tex", se star(* 0.10 ** 0.05 *** 0.01) drop(_cons) nomtitle f stats(N ymean year program, fmt( %9.0fc %7.3f %3s %3s) labels("Number of observations"  "Mean dependent variable" "Exam year FE" "Major FE")) cells(b(star pvalue(pboot) fmt(%9.3f)) seboot(par fmt(%9.3f))) replace collabels(none) nogaps booktabs label /// 
refcat(female " \\ \multicolumn{9}{l}{\textit{Panel A: Average annual wages (7-12 years after admission exam)}} \\", nolabel)

esttab l_max_annual_wage_712*  using "Output\Add_Table_AllP1_Annual_wages_`m'_ENEM_712_BS.tex", se star(* 0.10 ** 0.05 *** 0.01) drop(_cons) nomtitle f stats(N ymean year program , fmt( %9.0fc %7.3f %3s %3s) labels("Number of observations"  "Mean dependent variable" "Exam year FE" "Major FE")) cells(b(star pvalue(pboot) fmt(%9.3f)) seboot(par fmt(%9.3f))) append collabels(none) nogaps booktabs nonumber label ///
refcat(female " \\ \multicolumn{9}{l}{\textit{Panel B: Maximum annual wages (7-12 years after admission exam)}} \\", nolabel)


}

***** Erase temporary data sets
erase "Work Data/Temp_Bootstrap_AddTab.dta"
erase "Work Data/Temp_Wages_AddTable_BS.dta"






